function vtk3d_efield

fbase = 'example3d'
fnum  = 20

jx=48; jy=48; jz=16;   % box size


  fn = num2str(fnum,'%04d');

E = zeros(3,jx,jy,jz);

%-- electic field: Ex

  fname = [fbase, '.Ex.', fn];

  fid = fopen(fname, 'rb');
  v = fread(fid, jx*jy*jz, 'double');
  fclose(fid);

  v = reshape(v, jz,jy,jx);

  E(1,:,:,:) = flipxz(v);


%-- electic field: Ey

  fname = [fbase, '.Ey.', fn];

  fid = fopen(fname, 'rb');
  v = fread(fid, jx*jy*jz, 'double');
  fclose(fid);

  v = reshape(v, jz,jy,jx);

  E(2,:,:,:) = flipxz(v);


%-- electic field: Ez

  fname = [fbase, '.Ez.', fn];

  fid = fopen(fname, 'rb');
  v = fread(fid, jx*jy*jz, 'double');
  fclose(fid);

  v = reshape(v, jz,jy,jx);

  E(3,:,:,:) = flipxz(v);

whos E

  3*jx*jy*jz

%-- vtk header --
 
  vname = 'efield'
  fname = [fbase, '.efield.', fn, '.vtk'];

  fid=fopen(fname, 'wt');
  fprintf(fid, '# vtk DataFile Version 3.0\n');
  fprintf(fid, 'vtk output\n');
  fprintf(fid, 'BINARY\n');
  fprintf(fid, 'DATASET STRUCTURED_POINTS\n');
  fprintf(fid, 'DIMENSIONS %d %d %d\n', jx+1, jy+1, jz+1);
  fprintf(fid, 'SPACING 1 1 1\n');
  fprintf(fid, 'ORIGIN 0 0 0\n');
  fprintf(fid, 'POINT_DATA %d\n', jx*jy*jz);
  fprintf(fid, 'VECTORS %s float\n', vname);
  fprintf(fid, 'LOOKUP_TABLE default\n');
  fclose(fid);

%-- vtk data --

  fid=fopen(fname, 'ab', 'ieee-be');
  fwrite(fid, E, 'float');
  fclose(fid);



end


%---------------------------------------------------------

function fv = flipxz(v)

  [jz,jy,jx] = size(v);
  fv = zeros(jx,jy,jz);

  for j=1:jy
    a = v(:,j,:);
    a = reshape(a, jz, jx);
    fv(:,j,:) = a';
  end

end

%---------------------------------------------------------

